#if !defined(_Particle_Operator_test_)
#define _Particle_Operator_test_
#include <src/link/Qhull_Link.h>
#include <src/math/Operators/delta_particle.h>
#include <src/math/BoundaryCondition/BoundaryCondition.h>
#include <src/math/functor/GaussianFunctor.h>


template <class data, std::size_t Dim = data::Dim> 
struct particle_operator_fill_data  
{
	data operator()(std::size_t i)
	{
		OS_STATIC_CHECK(false); // you have to implement this test
	};
};

template <class data> 
struct particle_operator_fill_data<data,1>
{
	data operator()()
	{
		return 	InverseCumulativeFunctor<data>()(8);
	};
};

template <class data> 
struct particle_operator_fill_data<data,2> 
{
	data operator()() 
	{
		return 	InverseCumulativeFunctor<data>()(8,3);
	};
};


template <class particle_operator_type>
void particle_operator_example()
{
	typedef particle_operator_type::weight_set_type matrix_type;
	typedef particle_operator_type::argument_base_type argument_base_type;
	typedef particle_operator_type::result_type result_type;
	typedef particle_operator_type::state_variable_type state_variable_type;

	std::cout << std::endl << "----Start particle_operator_example---------------------" << std::endl ;
	std::cout << std::endl << " example for " << typeid(particle_operator_type).name()<< std::endl << std::endl ;

	argument_base_type* un = new argument_base_type();
	*un = particle_operator_fill_data<argument_base_type>()();
	std::cout << "A vector " << *un << std::endl;

	particle_operator_type operator1;
	operator1.set_state_variable(state_variable_type(un));

	std::cout << "Testing Algebra ..." << std::endl;

	particle_operator_type operator2;
	operator2 = operator1;
	OS_DYNAMIC_CHECK(operator2 == operator1, "particle_operator_example : test operator2 == operator1 failed");


	argument_base_type result = operator1(*un);
	std::cout << std::endl <<  typeid(particle_operator_type).name() << " applied to the vector " << result ;

	operator2 = operator1/2.;
	argument_base_type un_1 = operator2(*un);
	OS_DYNAMIC_CHECK(un_1 == result/2., "particle_operator_example : test un_1 == result/2. failed");
	operator2 = operator1+operator1;
	un_1 = operator2(*un);
	OS_DYNAMIC_CHECK(un_1 == result*2., "particle_operator_example : test un_1 == result*2. failed");
	operator2 = operator1-operator1;
	un_1 = operator2(*un);
	OS_DYNAMIC_CHECK(un_1 == 0.*result, "particle_operator_example : test un_1 == 0.*result failed");
	operator2 = 2.*(3.*operator1+operator1/2.);
	un_1 = operator2(*un);

	OS_DYNAMIC_CHECK(jmmath::norm2(un_1-7.*result)<THRESHOLD, "particle_operator_example : test un_1 == 7.*result failed");
	std::cout << "Output Operator 2.*(3.*operator1+operator1/2.) : " << std::endl<< operator2 << std::endl;
	particle_operator_type operator3;
	operator3 = operator2;


	OS_DYNAMIC_CHECK(operator3 == operator2, "particle_operator_example : test operator3 == operator2 failed");
	std::cout << std::endl << "----End particle_operator_example---------------------" << std::endl ;

};



void Particle_Operator_test();

#endif
